!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! Set up for the continuous Galerkin Stokes and Navier-Stokes codes.
!!
!! \n
!! 
!! This module collects all the initialization and finalization work
!! needed for the continuous Galerkin codes for the Stokes and
!! Navier-Stokes problems. This code would usually go in the main
!! program, but we collect it here to simplify the main and to share
!! it among various main programs.
!!
!! Regarding the input file, \c mod_cg_ns_setup_constructor reads the
!! <em>namelist</em> \c input from a file passed as argument. This
!! namelist collects all the input parameters which are common to all
!! the main programs. Each main program can then use an additional,
!! specialized input namelist which can be placed either in the same
!! file as the one used here or in a separate file. The \c input
!! <em>namelist</em> has the following structure:
!! <ul>
!!  <li> \c testname: the name of the chosen test case
!!  <li> \c basename: the base name for the output files, so that
!!  the output is written in files <em>basename-grid.octxt</em>,
!!  <em>basename-base.octxt</em>,
!!  <em>basename-results.octxt</em> and so on
!!  <li> \c ku: degree of the velocity finite element space
!!  <li> \c kb: degree of the bubble function(s) for the velocity
!!  finite element space, all the bubble functions with index ranging
!!  from <tt>kb(1)</tt> to <tt>kb(2)</tt> are included, use negative
!!  values to not include any bubble degree of freedom
!!  <li> \c kp: degree of the pressure finite element space
!!  <li> \c pressure_stabilization: set this to <tt>.true.</tt> to
!!  include a pressure stabilization allowing equal order
!!  approximation for velocity and pressure
!!  <li> \c p_stab_type: type of the pressure stabilization
!!  term (only used if \c pressure_stabilization is <tt>.true.</tt>);
!!  see \c mod_cg_ns_linsystem for the available stabilizations.
!!  <li> \c p_stab_coeff: coefficient of the pressure stabilization
!!  term (only used if \c pressure_stabilization is <tt>.true.</tt>)
!!  <li> \c write_grid: set this to <tt>.true.</tt> to write the grid
!!  to a file
!!  <li> \c grid_file: path to the grid file (relative to the working
!!  directory)
!!  <li> \c nbound: number of boundary regions (if, for instance,
!!  there are three boundary regions with indices 1, 2, 4, one should
!!  set <tt>nbound=4</tt>; the fact that the index 3 is never used is
!!  not a problem)
!!  <li> \c cbc_type: this is an <tt>2*nbound</tt> integer array where
!!  the first \c nbound values are 1 for Dirichlet bcs and 2 for
!!  Neumann bcs, and the last nbound values specify the kind of
!!  boundary condition as discussed in \c mod_cgstokes_locmat.
!!  <li> \c add_div_lagmultiplier: set this to <tt>.true.</tt> to
!!  include an additional scalar unknown
!!  \f{displaymath}{
!!   \lambda = \int_{\Omega} \nabla\cdot\underline{u} \, dx,
!!  \f}
!!  which is the degree of freedom associated with the Lagrange
!!  multiplier of the constraint
!!  \f{displaymath}{
!!   \int_{\Omega} p \, dx = 0;
!!  \f}
!!  this option is necessary when all the boundary conditions are of
!!  Dirichlet type, since otherwise the pressure can not be defined
!!  uniquely (this additional variable also eliminates the constraint
!!  of zero total flux on the Dirichlet datum)
!!  <li> \c linsolver: linear solver, currently implemented are
!!  "umfpack", "mumps" and "iterative" (see then \c mod_linsolver and
!!  \c mod_iterativesolvers for the available iterative solvers). When
!!  using the umfpack solver, the following variables are also
!!  relevant:
!!  <ul>
!!   <li> \c umfpack_print_level: this value is used to set
!!   <tt>umf_contr(umfpack_prl)</tt>, which determines the verbosity
!!   of the umfpack subroutines as discussed in the relative
!!   documentation.
!!  </ul>
!!  When using the iterative solver, the following variables are also
!!  relevant:
!!  <ul>
!!   <li> \c itsolver: the solver to be used
!!   <li> \c itsol_abstol: set this to <tt>.true.</tt> to use an
!!   absolute tolerance in the iterative solver (default is relative)
!!   <li> \c itsol_toll: tolerance for the iterative solver
!!   <li> \c preconditioner: specifies the preconditioner (see \c
!!   mod_lhs for all the available options)
!!   <li> \c precon_ndiags: if
!!   <tt>preconditioner.eq.'multidiagonal'</tt>, this value specifies
!!   how many diagonals of the complete matrix form the
!!   preconditioner; a negative value means using the whole matrix
!!  </ul>
!!  <li> \c write_sys: set this to <tt>.true.</tt> to include the
!!  linear system in the output file (the file of course will be
!!  significantly bigger)
!!  <li> \c compute_error_norms: set this to <tt>.true.</tt> to
!!  compute the error norms; this assumes that the exact solution is
!!  available in the corresponding test case file
!!  <li> \c err_extra_deg: when computing error norms use quadrature
!!  rules of degrees \f$default\,\,degree + err_{extra\,\,deg}\f$,
!!  where \f$default\,\,degree\f$ is the degree used in the
!!  computation of the finite element system.
!! </ul>
!<----------------------------------------------------------------------
module mod_cg_ns_setup

!-----------------------------------------------------------------------

 use mod_utils, only: &
   t_realtime, my_second

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_constructor, &
   mod_kinds_destructor,  &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_constructor, &
   mod_fu_manager_destructor,  &
   new_file_unit

 use mod_octave_io, only: &
   mod_octave_io_constructor, &
   mod_octave_io_destructor

 use mod_sympoly, only: &
   mod_sympoly_constructor, &
   mod_sympoly_destructor,  &
   nll_sympoly

 use mod_octave_io_sympoly, only: &
   mod_octave_io_sympoly_constructor, &
   mod_octave_io_sympoly_destructor

 use mod_linal, only: &
   mod_linal_constructor, &
   mod_linal_destructor

 use mod_perms, only: &
   mod_perms_constructor, &
   mod_perms_destructor

 use mod_octave_io_perms, only: &
   mod_octave_io_perms_constructor, &
   mod_octave_io_perms_destructor

 use mod_sparse, only: &
   mod_sparse_constructor, &
   mod_sparse_destructor

 use mod_octave_io_sparse, only: &
   mod_octave_io_sparse_constructor, &
   mod_octave_io_sparse_destructor

 use mod_mpi_utils, only: &
   mod_mpi_utils_constructor, &
   mod_mpi_utils_destructor,  &
   mpi_comm_world, &
   mpi_init, mpi_finalize, &
   mpi_comm_size, mpi_comm_rank

 use mod_output_control, only: &
   mod_output_control_constructor, &
   mod_output_control_destructor,  &
   elapsed_format, &
   base_name

 use mod_master_el, only: &
   mod_master_el_constructor, &
   mod_master_el_destructor,  &
   t_lagnodes, lag_nodes, &
   add_nodes

 use mod_grid, only: &
   mod_grid_constructor, &
   mod_grid_destructor,  &
   t_grid, t_ddc_grid,   &
   new_grid, clear,      &
   write_octave

 use mod_bcs, only: &
   mod_bcs_constructor, &
   mod_bcs_destructor,  &
   t_bcs, new_bcs, clear, &
   write_octave

 use mod_numquad, only: &
   mod_numquad_constructor, &
   mod_numquad_destructor

 use mod_state_vars, only: &
   mod_state_vars_constructor, &
   mod_state_vars_destructor,  &
   c_stv

 use mod_linsolver, only: &
   mod_linsolver_constructor, &
   mod_linsolver_destructor,  &
   c_linpb, c_itpb, c_mumpspb

 use mod_base, only: &
   mod_base_constructor, &
   mod_base_destructor,  &
   t_base, clear, &
   write_octave

 use mod_testcases, only: &
   mod_testcases_constructor, &
   mod_testcases_destructor

 use mod_fe_spaces, only: &
   mod_fe_spaces_constructor, &
   mod_fe_spaces_destructor,  &
   cg ! scalar f.e. space

 use mod_cgdofs, only: &
   mod_cgdofs_constructor, &
   mod_cgdofs_destructor,  &
   t_cgdofs, t_ddc_cgdofs, &
   new_cgdofs,           &
   clear,                &
   write_octave

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_cg_ns_setup_constructor, &
   mod_cg_ns_setup_destructor,  &
   mod_cg_ns_setup_initialized, &
   ! miscellaneous
   max_char_len,   &
   ! Grid
   grid, ddc_grid, bcs, &
   ! FE space
   pressure_stabilization, add_div_lagmultiplier, &
   ku, kp,                      &
   p_stab_type, p_stab_coeff,   &
   uref_nodes, pref_nodes, &
   ubase,     pbase,     &
   udofs,     pdofs,     &
   ddc_udofs, ddc_pdofs, &
   ! Linear system
   mpi_nd, mpi_id, linsolver, &
   t_nodal_field, t_linpb_mumps, linpb, &
   ! error computation
   compute_error_norms, &
   err_extra_deg,  &
   ! IO
   write_sys,      &
   out_file_nml_name

 private
 save

!-----------------------------------------------------------------------

 !----------------------------------------------------------------------
 ! Public members

 ! miscellaneous
 integer, parameter :: max_char_len = 1000
 ! FE space
 logical, protected :: pressure_stabilization, add_div_lagmultiplier
 integer, protected :: ku, kp, kb(2)
 character(len=max_char_len), protected :: p_stab_type
 real(wp), protected :: p_stab_coeff
 type(t_lagnodes), allocatable, protected :: uref_nodes(:), pref_nodes(:)
 type(t_base), protected :: ubase, pbase
 type(t_cgdofs    ), protected :: udofs, pdofs
 type(t_ddc_cgdofs), protected :: ddc_udofs, ddc_pdofs
 ! Grid
 type(t_grid    ), target :: grid
 type(t_ddc_grid)         :: ddc_grid
 type(t_bcs), target :: bcs
 ! Linear system
 integer :: mpi_nd, mpi_id
 class(c_linpb), allocatable :: linpb
 !> Nodal field used in the solution of the linear system
 type, extends(c_stv) :: t_nodal_field
  real(wp), allocatable :: x(:) ! nodal values
 contains
  procedure, pass(x) :: incr => nodal_field_incr
  procedure, pass(x) :: tims => nodal_field_tims
  procedure, pass(z) :: copy => nodal_field_copy
  !!> The following procedures are introduced for efficiency reasons
  !procedure, pass(x) :: inlt => nodal_field_inlt
  procedure, pass(x) :: source => nodal_field_source
  procedure, pass(x) :: source_vect => nodal_field_source_vect
 end type t_nodal_field
 !> Linear problem (MUMPS viewpoint)
 type, extends(c_mumpspb) :: t_linpb_mumps
 contains
  procedure, pass(s) :: xassign => mumps_xassign
 end type t_linpb_mumps
 ! Error computation
 logical, protected :: compute_error_norms
 integer, protected :: err_extra_deg
 ! IO
 logical, protected :: write_sys
 character(len=*), parameter :: out_file_nml_suff = '-nml.octxt'
 character(len=max_char_len+len(out_file_nml_suff)), protected :: &
   out_file_nml_name
 ! local
 logical, protected :: mod_cg_ns_setup_initialized = .false.
 !----------------------------------------------------------------------

 !----------------------------------------------------------------------
 ! Private members

 ! grid
 integer :: nbound
 character(len=max_char_len) :: grid_file, cbc_type
 ! test case
 character(len=max_char_len) :: testname
 ! linear system solver
 logical :: itsol_abstol
 integer :: umfpack_print_level, precon_ndiags
 real(wp) :: itsol_toll
 character(len=max_char_len) :: linsolver, itsolver, &
   preconditioner
 ! IO
 logical :: write_grid
 character(len=max_char_len) :: basename
 ! local
 character(len=*), parameter :: this_mod_name = 'mod_cg_ns_setup'
 !----------------------------------------------------------------------

 ! Input namelist
 namelist /input/ &
   ! test case
   testname, &
   ! output base name
   basename, &
   ! base
   ku, kb, kp,  &
   pressure_stabilization, p_stab_type, p_stab_coeff, &
   ! grid
   write_grid,  &
   grid_file,   &
   ! boundary conditions
   nbound, cbc_type, &
   add_div_lagmultiplier, &
   ! linear system solver
   linsolver,   &
   ! UMFPACK control variables (used only if linsolver.eq."umfpack")
   umfpack_print_level, &
   ! It. sol. control variables (used only if linsolver.eq."iterative")
   itsolver, itsol_abstol, itsol_toll, &
   !  precon_ndiags can be larger than the size of the matrix, in
   !  which case the whole matrix is used as preconditioner
   preconditioner, precon_ndiags, &
   ! results
   write_sys,   &
   ! errors
   compute_error_norms, err_extra_deg

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 !> \details This constructor, besides doing various initializations,
 !! reads the input file. The input file is taken from the command
 !! line arguments, if present, otherwise \c input_file_name_def is
 !! used. The optional argument \c input_file_name_act can be used to
 !! return the name of the input file chosen after checking the
 !! command line parameters.
 !<
 subroutine mod_cg_ns_setup_constructor( input_file_name_def, &
                                         input_file_name_act )
  character(len=*), intent(in) :: input_file_name_def
  character(len=*), intent(out), optional :: input_file_name_act
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

  ! grid
  integer, allocatable :: bc_type_t(:,:)
  ! FE space
  integer :: nbubs ! number of bubbles
  ! IO variables
  character(len=max_char_len) :: input_file_name
  character(len=*), parameter :: out_file_grid_suff = '-grid.octxt'
  character(len=10000+len(out_file_grid_suff)) :: out_file_grid_name
  character(len=*), parameter :: out_file_base_suff = '-base.octxt'
  character(len=10000+len(out_file_base_suff)) :: out_file_base_name
  ! Auxiliary variables
  integer :: fu, ierr
  character(len=10*max_char_len) :: message(1)
  ! timing
  real(t_realtime) :: t0, t1

   !Consistency checks ---------------------------
   if(mod_messages_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_cg_ns_setup_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   !---------------------------------------------------------------------
   ! Setup of all the MPI processes
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! Preliminary initializations
   call mod_kinds_constructor()

   call mpi_init(ierr)
   call mod_mpi_utils_constructor()
   call mpi_comm_size(mpi_comm_world,mpi_nd,ierr)
   call mpi_comm_rank(mpi_comm_world,mpi_id,ierr)

   call mod_fu_manager_constructor()

   call mod_octave_io_constructor()

   call mod_sympoly_constructor()
   call mod_octave_io_sympoly_constructor()

   call mod_linal_constructor()

   call mod_perms_constructor()
   call mod_octave_io_perms_constructor()

   call mod_sparse_constructor()
   call mod_octave_io_sparse_constructor()
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! Read input file
   if(command_argument_count().gt.0) then
     call get_command_argument(1,value=input_file_name)
   else
     input_file_name = input_file_name_def
   endif

   call new_file_unit(fu,ierr)
   open(fu,file=trim(input_file_name), &
        status='old',action='read',form='formatted',iostat=ierr)
   read(fu,input)
   close(fu,iostat=ierr)
   ! If there are more processes, each of them needs its own basename
   ! and grid
   if(mpi_nd.gt.1) then
     write(basename, '(a,a,i3.3)') trim(basename),'-P',mpi_id
     write(grid_file,'(a,a,i3.3)') trim(grid_file),'.',mpi_id
   endif
   ! echo the input namelist
   out_file_nml_name = trim(basename)//out_file_nml_suff
   open(fu,file=trim(out_file_nml_name), &
        status='replace',action='write',form='formatted',iostat=ierr)
   write(fu,input)
   close(fu,iostat=ierr)

   allocate(bc_type_t(nbound,2))
   read(cbc_type,*) bc_type_t ! transposed, for simplicity

   if(present(input_file_name_act)) &
     input_file_name_act = trim(input_file_name)
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! output control setup
   call mod_output_control_constructor(basename)
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! linear solver (the solver will be defined once we know the matrix)
   call mod_state_vars_constructor()
   call mod_linsolver_constructor()
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! Construct the grid (without the continuous dofs, which will be
   ! computed later once the FE basis has been defined)
   t0 = my_second()

   call mod_master_el_constructor()

   call mod_grid_constructor()
   call new_grid(grid,trim(grid_file) , ddc_grid,mpi_comm_world)

   ! write the octave output
   if(write_grid) then
     out_file_grid_name = trim(base_name)//out_file_grid_suff
     call new_file_unit(fu,ierr)
     open(fu,file=trim(out_file_grid_name), &
          status='replace',action='write',form='formatted',iostat=ierr)
      call write_octave(grid    ,'grid'    ,fu)
      call write_octave(ddc_grid,'ddc_grid',fu)
     close(fu,iostat=ierr)
   endif

   t1 = my_second()
   write(message(1),elapsed_format) &
     'Completed grid construction: elapsed time ',t1-t0,'s.'
   call info(this_sub_name,this_mod_name,message(1))
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! Collect the information concerning the boundary conditions
   call mod_bcs_constructor()
   call new_bcs(bcs,grid,transpose(bc_type_t) , ddc_grid,mpi_comm_world)
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! Define the finite element space
   t0 = my_second()

   call mod_numquad_constructor()

   call mod_base_constructor()

   call mod_fe_spaces_constructor()

   call mod_cgdofs_constructor()

   ! nodal Lagrangian basis
   call lag_nodes(uref_nodes,grid%me%d,ku)
   if(kb(2).ge.0) then ! add the bubbles
     if(kb(1).lt. 0   ) kb(1) = 0
     if(kb(1).gt.kb(2)) call error(this_sub_name,this_mod_name, &
                  'The value kb(1) cannot be larger than kb(2).')
     ubase = cg(grid%me,ku,uref_nodes,b_k=kb)
     ! The easiest way to deal with the bubble functions is now by
     ! adding one internal node for each bubble function. However, this
     ! is only a useful hack, since the bubble functions are not
     ! Lagrangian and they are not associated with any node. The node
     ! that we associate with them are internal and conventionally
     ! placed at the barycenter.
     ! First we need to determine the number of bubbles: see cg
     nbubs = nll_sympoly(grid%d,kb(2))
     if(kb(1).gt.0) nbubs = nbubs - nll_sympoly(grid%d,kb(1)-1)
     call add_nodes( uref_nodes(grid%d), & ! add internal nodes
            spread( & ! new nodes are in the barycenter
            spread(1.0_wp/(grid%d+1),1,grid%d),2,nbubs ) &
                   )
   else
     ubase = cg(grid%me,ku,uref_nodes)
   endif
   call new_cgdofs( udofs,uref_nodes,grid,bcs, &
             mpi_comm_world,ddc_grid,ddc_udofs )

   call lag_nodes(pref_nodes,grid%me%d,kp)
   pbase = cg( grid%me , kp , pref_nodes ,                 &
               deg=ubase%deg , xig=ubase%xig , wg=ubase%wg )
   call new_cgdofs( pdofs,pref_nodes,grid,bcs, &
             mpi_comm_world,ddc_grid,ddc_pdofs )

   ! write the octave output
   out_file_base_name = trim(base_name)//out_file_base_suff
   call new_file_unit(fu,ierr)
   open(fu,file=trim(out_file_base_name), &
        status='replace',action='write',form='formatted',iostat=ierr)
     call write_octave(ubase,'ubase',fu)
     call write_octave(pbase,'pbase',fu)
     call write_octave(uref_nodes,'uref_nodes',fu)
     call write_octave(pref_nodes,'pref_nodes',fu)
     if(write_grid) then
       call write_octave(bcs  ,'bcs'  ,fu)
       call write_octave(udofs,'udofs',fu)
       call write_octave(pdofs,'pdofs',fu)
       call write_octave(ddc_udofs,'ddc_udofs',fu)
       call write_octave(ddc_pdofs,'ddc_pdofs',fu)
     endif
   close(fu,iostat=ierr)

   t1 = my_second()
   write(message(1),elapsed_format) &
     'Completed FE basis construction: elapsed time ',t1-t0,'s.'
   call info(this_sub_name,this_mod_name,message(1))
   !---------------------------------------------------------------------

   !---------------------------------------------------------------------
   ! Test case setup
   call mod_testcases_constructor(testname,grid%me%d)
   !---------------------------------------------------------------------

   deallocate(bc_type_t)
   mod_cg_ns_setup_initialized = .true.
 end subroutine mod_cg_ns_setup_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_cg_ns_setup_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'

  integer :: ierr
   
   !Consistency checks ---------------------------
   if(mod_cg_ns_setup_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call mod_testcases_destructor()

   call clear( bcs )
   call mod_bcs_destructor()

   call clear( ubase );     call clear( pbase )
   call clear( udofs );     call clear( pdofs )
   call clear( ddc_udofs ); call clear( ddc_pdofs )
   deallocate( uref_nodes , pref_nodes )

   call mod_fe_spaces_destructor()

   call mod_base_destructor()
  
   call mod_numquad_destructor()

   call mod_cgdofs_destructor()
  
   call clear( grid )
   call clear( ddc_grid )
   call mod_grid_destructor()
  
   call mod_master_el_destructor()
  
   call mod_linsolver_destructor()
   call mod_state_vars_destructor()

   call mod_output_control_destructor()

   call mod_octave_io_sparse_destructor()
   call mod_sparse_destructor()
  
   call mod_octave_io_perms_destructor()
   call mod_perms_destructor()
  
   call mod_linal_destructor()
  
   call mod_octave_io_sympoly_destructor()
   call mod_sympoly_destructor()
  
   call mod_octave_io_destructor()
  
   call mod_fu_manager_destructor()
  
   call mod_mpi_utils_destructor()
   call mpi_finalize(ierr)

   call mod_kinds_destructor()
   !--------------------------------------------------------------------

   mod_cg_ns_setup_initialized = .false.
 end subroutine mod_cg_ns_setup_destructor

!-----------------------------------------------------------------------

 subroutine nodal_field_incr(x,y)
  class(c_stv),         intent(in)    :: y
  class(t_nodal_field), intent(inout) :: x

   select type(y); type is(t_nodal_field)
    x%x = x%x + y%x
   end select
 end subroutine nodal_field_incr

!-----------------------------------------------------------------------

 subroutine nodal_field_tims(x,r)
  real(wp),             intent(in)    :: r
  class(t_nodal_field), intent(inout) :: x

   x%x = r*x%x
 end subroutine nodal_field_tims

!-----------------------------------------------------------------------

 subroutine nodal_field_copy(z,x)
  class(c_stv),         intent(in)    :: x
  class(t_nodal_field), intent(inout) :: z

   select type(x); type is(t_nodal_field)
    z%x = x%x
   end select
 end subroutine nodal_field_copy

!-----------------------------------------------------------------------

 !> Workaround for ifort compiler bug (see \c mod_state_vars)
 subroutine nodal_field_source(y,x)
  class(t_nodal_field), intent(in)  :: x
  class(c_stv), allocatable, intent(out) :: y

   allocate(t_nodal_field::y)
   select type(y); type is(t_nodal_field)
    allocate(y%x(size(x%x)))
    y%x = x%x
   end select
 end subroutine nodal_field_source
 subroutine nodal_field_source_vect(y,x,m)
  integer, intent(in) :: m
  class(t_nodal_field), intent(in)  :: x
  class(c_stv), allocatable, intent(out) :: y(:)

  integer :: i

   allocate(t_nodal_field::y(m))
   select type(y); type is(t_nodal_field)
    do i=1,m
      allocate(y(i)%x(size(x%x)))
      y(i)%x = x%x
    enddo
   end select
 end subroutine nodal_field_source_vect

!-----------------------------------------------------------------------

 subroutine mumps_xassign(x,s,mumps_x)
  real(wp),             intent(in)    :: mumps_x(:)
  class(t_linpb_mumps), intent(inout) :: s
  class(c_stv),         intent(inout) :: x

   select type(x); type is(t_nodal_field)
    x%x = mumps_x
   end select
 end subroutine mumps_xassign

!-----------------------------------------------------------------------

end module mod_cg_ns_setup

